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We show that the Hamiltonian mean field (HMF) model describes the equilibrium behavior 
of a system of long pendula with flat bobs that are coupled through long-range interactions 
(charged or self gravitating). We solve for the canonical partition function in the coordinate 
frame of the pendula angles. The Hamiltonian in the angles coordinate frame looks similar 
to the form of the HMF model but with the inclusion of an index dependent phase in the 
interaction term. We also show interesting non-equilibrium behavior of the pendula angles, 
namely that a quasistationary clustered state can exist when pendula angles are initially 


ordered by their index. 

I. INTRODUCTION 

Systems with long-range interactions are a 
source of unique problems in the field of statisti¬ 
cal mechanics and thermodynamics. This is due 
to several properties of long-range systems which 
fall outside of the conditions normally needing 
to be satisfied when applying the methodolo¬ 
gies of thermodynamics. Simply from the words 
“long-range” the first infringement can be de¬ 
duced, that long-range systems are not addi¬ 
tive. If two systems with short-range interac¬ 
tions are brought together to form a larger sys¬ 
tem then the energy difference between the con¬ 
glomerate system and the sum of its constituents 
is the new potential energy from the boundary 
between them. In the thermodynamic limit, the 
potential energy of the boundary is small com¬ 
pared to the bulk and can be neglected, mak¬ 
ing short-range systems additive. In the case of 
long-range interactions, one particle will feel a 
significant potential created by every other par¬ 
ticle, so the additional potential energy of two 
systems added together does not scale as the 
boundary but in a more complicated way that 
depends on the specific nature of the interactions 
jQj]. Directly related to the lack of additivity is 
the fact that systems with long-range interac¬ 
tions are not extensive because their energy di¬ 
verges in the thermodynamic limit |2]. Although 

* oweenm@gmail.com 


these characteristics compel cautious use of the 
usual tools of statistical mechanics, they are also 
the source of many interesting dynamical and 
statistical features. Depending on the system of 
interest, such features include canonical and mi- 
crocanonical ensemble inequivalence and related 
negative specific heat 0, quasistationary states 
(different than metastable states which lie on 
local extrema of equilibrium potentials) whose 
lifetimes increase with the number of particles 
@3, an interesting dependence of the largest Lya¬ 
punov exponent on particle number in a long- 
range Fermi-Pasta-Ulam model (5], and sponta¬ 
neous creation of macroscopic structures in non¬ 
equilibrium states [6|. In some cases, long-range 
interactions can greatly simplify problems. For 
instance, mean field models depend on one of two 
premises: (i) interactions are short-range but the 
system is embedded in a space of infinite dimen¬ 
sion so that all bodies in the system are nearest 
neighbors, or (ii) interactions are infinitely long. 

For some time, the primary motivation 
for the study of long-range interactions was 
to understand galaxies, galaxy clusters and 
the general thermodynamic properties of self- 
gravitating systems. Aside from mean field mod¬ 
els, interest has further built since the obser¬ 
vation of modified scattering lengths in Bose- 
Einstein Condensates (BEC) through the use of 
Feshbach resonances |7j. Using this technique, a 
BEC can be made to be almost non-interacting 
by tuning the scattering length to zero. One 
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could even tune the scattering length to a neg¬ 
ative value, making the BEC collapse. More re¬ 
cently, O’dell et al. [8] has shown that it may be 
possible to produce an attractive 1/r potential 
between atoms in a BEC by applying an “ex¬ 
tremely off resonant” electromagnetic field. This 
has opened the possibility of creating table-top 
methods which physically model aspects of cos¬ 
mological behavior on a laboratory scale, as well 
as the possible development of entirely new dy¬ 
namics in BEC. 

The challenges in understanding long-range 
systems drive the development of solvable mod¬ 
els that could help better explain some of the 
aforementioned phenomena. Carnpa et al. [9] 
have recently published a collection of impor¬ 
tant solvable models. One particularly signifi¬ 
cant model, which is important to this work, is 
the Hamiltonian Mean Field (HMF) XY spin 
model |6], often written in the form 


H = E| + 5V £ [i-<*»(*-«;)]. (i) 

i=1 i,j =1 

where 9{ is the angular position of the i th parti¬ 
cle (spin), as shown in Fig. [lj and pi is its con¬ 
jugate angular momentum. The HMF model is 
generally used to describe two different classes 
of systems: 1) a mean field XY classical spin 
model, and 2) a one dimensional periodic system 
of itinerant particles with long-range interac¬ 
tions. Though the connection between the HMF 
model and the second class of systems mentioned 
could be thought of as contrived given the sim¬ 
plifications under which the model is realized, it 
has been shown that the model produces useful 
insights into how non-neutral plasmas and self 
gravitating systems behave [6]. 

In this paper, we study the dynamics of an 
array of N pendula with long-range interacting 
bobs. By considering long pendula with flat bobs 
undergoing small oscillations and having parallel 
planes of rotation, we produce a model related 
to the HMF model through a coordinate trans¬ 
formation. The transformation introduces a de¬ 
pendence on the indices of the particle labels. A 
cartoon of the physical picture is shown in Fig. [lj 


The index dependence in the Hamiltonian, that 



Figure 1. N pendulum system with parallel planes 
of rotation. The i th pendulum angle at some time t 
is 9i{t). 

will be described in detail in the next section, is a 
consequence of the pendula pivots being slightly 
offset from one another and appears as a phase 
in the cosine term of the HMF model. It inspires 
the investigation of non-equilibrium “repulsive” 
behavior in the angle coordinate frame where we 
find an interesting quasistationary state when 
the angles of the pendula are initially ordered 
according to their indices. We find the clustered 
positions in the usual HMF coordinate frame (bi¬ 
clusters) , but in the angle coordinate frame clus¬ 
tering is only found for the initially ordered an¬ 
gles and, unlike the biclusters, these are clearly 
quasistationary states. A quasistationary state 
is defined as a dynamical state that can persist 
for a length of time which goes to infinity as 
the thermodynamic limit is approached [9j. In 
addition to discussing the clustered angle states 
exhibited by the system, we also solve for the 
canonical partition function in the pendulum an¬ 
gle coordinate frame, finding that in equilibrium 
with a heat bath, the probability distributions of 
the angles can be described by the original HMF 
model. This finding is similar to the work done 
by HD| on a model sometimes called the HMF a- 
model. In the HMF a model, a 1/V“- dependence 
between the classical spins is introduced [a tm- 
llij . where r l3 is the distance between the and 
j th spins on a lattice. Though the physical mo¬ 
tivations behind studying these various models 
can be very different, it is interesting that their 
equilibrium behavior is the same or nearly the 
same. We believe that the work in this paper 
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further suggests that the HMF model universally 
describes an entire class of long-range interact¬ 
ing systems in equilibrium. 


terms of c pi, the positions can be rewritten as 


Xi = 


27 vi 

W 


+ Pi- 


(4) 


B. Density Approximation 

II. THE MODEL 


A. Coordinates 


In Fig. [lj we show an array of pendula ro¬ 
tating in the same plane with bobs that inter¬ 
act through a long-range potential. If we con¬ 
sider the case where all the pendula only undergo 
small oscillations, we may write the horizontal 
location of the z th particle, Xi, as x t = id + £9i , 
where d is the distance between the pivots of 
neighboring pendula and £ is the length of each 
pendulum. The small 0 regime makes the prob¬ 
lem one dimensional in x. We choose periodic 
boundary conditions and rescale the system by 
2tt/N d so that 

2vr 

1 -> Wi x (2) 

making the total system length a dimensionless 
27 t where N is the number of particles in one 
period. We will refer to a periodic space with 
length 2ir as a unit circle. The position of the 
i th particle (bob) is now 


Xi = 


2ni 2ir l 

aT + n~/' 


(3) 


We have not yet explicitly stated the phys¬ 
ical mechanism through which the bobs inter¬ 
act. Connecting the interactions with specific 
physical motivations should be discussed with 
some discretion because the development of the 
model leaves these motivations up to some free¬ 
dom of interpretation. Imagine that the bobs all 
carry some charge. We will not distinguish be¬ 
tween particles in any other way than their in¬ 
dices, so in the case where all particles carry the 
same charge, repulsive behavior is expected. On 
the other hand one could make the bobs attract 
one another, which could be thought of as the 
self-gravitating case. To be solvable, the model 
requires some simplifications. For the sake of 
brevity we will speak of the particle charge or 
mass density as the “density”. 

The approximation that we invoke is simi¬ 
lar to that used when justifying the HMF model 
(Eq. Q) to describe free particles in a one¬ 
dimensional ring mm- The distribution of the 
bobs is such the mass density, p(x), is given by 

N i 

P( x ) =^2 6 ( x ~ x i) ~ ( 5 ) 

2=1 


For reasonable choices of £ and d {i/d « N), 
the second term on the RHS is suitably small 
such that the Hamiltonian can be written with 
terms that are quadratic in 9. However, we are 
primarily interested in a regime where i/d —> 
oo as the thermodynamic limit is approached. 
Physically this corresponds to the small oscil¬ 
lations of very long pendula with suspension 
points that are close together compared to their 
lengths. In order to simplify the calculations 
that follow, we define (pi to be the last term 
on the RHS of Eq. Q, namely (pi = 2nl6i/Nd. 
Given the choice of large i/d, pi can take any 
value in the range [0, 2ir). This is only true be¬ 
cause i/d is large, not because the 9iS are. In 


The constant \/2tx subtracted from the delta 
function is necessary to produce a meaningful ex¬ 
pression for the potential and corresponds to 
the inclusion of a neutralizing (of opposite sign) 
homogeneous background density. Restricting 
the problem further to that of solving Poisson’s 
equation for a one-dimensional potential phys¬ 
ically amounts to choosing large and flat bob 
geometries oriented with their smallest axis par¬ 
allel to the x axis. Writing the delta function 
as a cosine Fourier series, Poisson’s equation be¬ 
comes 

N oo 

V 2 $(x) = — y cos [n(x — Xi)]. (6) 

7r z J ' 

2=1 72=1 
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The parameter 7 contains the particle (bob) 
charge or mass and becomes the interaction 
strength in the Hamiltonian. We can see that 
the zeroth-order term in the Fourier series can¬ 
celed the constant neutralizing background that 
was superficially added. 

The most important simplification in this pa¬ 
per is truncating the sum of the Fourier coeffi¬ 
cients used to represent the delta function after 
the n = 1 coefficient. Antoni et al. defend the 
truncation by asserting that the “large scale col¬ 
lective properties” do not greatly change when 
higher order terms of the sum (including interac¬ 
tions at the smaller length scales) are included, 
and discuss the consequences of the approxima¬ 
tion in some detail [ 6 ]. The simplification also 
warrants a brief discussion of the way that it 
could be physically interpreted. The truncation 
of the sum is equivalent to smearing out the den¬ 
sity of each particle over the system so that it is 
peaked at its given location, Xi, but also having 
a negative density peak on the opposite side of 
the unit circle. This could be thought of as dou¬ 
bling the number of particles and enforcing that 
each particle has a negative partner that always 
remains on the opposing side of the unit circle. 
After this doubling, the now nebulous masses are 
dispersed such that a pair’s density is described 
by a cosine function with the positive peak cen¬ 
tered at Xi. 


ble, respectively. The force that the j th particle 
experiences when cj)j and all fa are zero is given 
by 


— V<F (xj) 



2 ?r (j - i ) 
N 



(8) 

The sum JT sin [2n(j — i)/N] equals zero for any 
j, so ci must be zero. Integrating once more to 
obtain the potential yields 



To determine C 2 we stipulate that if all fa = 0, 
then 41(0) = 0. Inserting Eq. Q (or Eq. Q) for 
Xi yields 



( 10 ) 


The sum over the cosine is zero, therefore C 2 = 0 
and we can now write the potential energy of the 
j th particle as 




7 

7r 


N 


^2 cos 
2=1 


2 ?r(j - i) 
N 


4*1 ■ 

(ID 


C. Solving Poisson’s Equation 


D. The Hamiltonian 


Integrating Poisson’s equation once, we ob¬ 
tain: 


The Hamiltonian can be written as 


N 

V<f>(x) = — {sin (x — Xi) + ci} . 
7r z ' 

2=1 


(7) 


H = H 0 + H I , 

where Hq is the kinetic energy piece 


In order to determine the constant ci from the 
integration, the physical picture should be exam¬ 
ined. A sensible requirement is that when all of 
the bobs are hanging at their equilibrium posi¬ 
tions, directly below their pivot (all fa = 9i = 0 ), an< l 
the net force experienced by any bob is zero. 

This is a valid requirement if the bobs are at- } 
tractive or repulsive, the only difference being 
that the configuration would be unstable or sta- 


N « 2 


2=1 


N 

y 

2N ^ 
1,3 =1 


cos 


27T ji-j) 
N 


( 12 ) 


(13) 


4>j 

( 14 ) 
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is the interaction piece, so 


N 


* = E 


i= 1 


Pj_ 

2 


7 

2 N 


N 

E 


cos 


2tt(z - j) 
N 


+ 4>i ~ ' 


(15) 

The mass of the bobs has been set to unity, 7 is 
the interaction strength, a factor of 1/2 accounts 
for the double counting, and the 1 / 7 r coefficient 
in the potential energy has been absorbed into 
7. The factor of 1/N is a rescaling of the poten¬ 
tial energy that ensures that as the thermody¬ 
namic limit is approached, the potential energy 
of the system does not diverge. The 1/N scal¬ 
ing is known as the Kac prescription [16] . The 
Kac serves to keep both the energy and entropy 
of a system proportional to the number of parti¬ 
cles in the system, an important prerequisite for 
phase transitions [5]. 


E. Relationship to the HMF model and the 
Spin Interpretation 

Due to the simple bijective relationship be¬ 
tween Xi and fa one can simply solve the equa¬ 
tions of motion for the HMF model and find 
the dynamics for (p t via the coordinate transform 
Xi —> fa. Previously it was mentioned that the 
HMF model is used to describe free particles on 
a ring with long-range repulsion or attraction, 
as well as describing a classical XY spin model. 
The Oi played the role of either the position of 
the i th particle on the ring or the orientation of 
the i th spin. Therefore, it is interesting to specu¬ 
late about the type of spin system the model de¬ 
scribes in the (pi picture. Thus far, the rescaled 
angle fa = 2ttIQ 1 /N(1 describes the distance of 
a pendulum bob from the point directly below 
its pivot, but it could also be interpreted as the 
orientation of spin. In the spin interpretation 
of Eq. (15), the potential energy of the i th and 


librium with a heat bath by solving the parti¬ 
tion function in the fa coordinate frame. In the 
process of simplifying the Hamiltonian to solve 
for the partition function, we will find expres¬ 
sions of the form cos fa and sin 4>i which we talk 
about as the horizontal and vertical components 
of a magnetization rhi = (cos<fo,sin<fo). It could 
easily be stated that in the spin analogy, the 
<f>i are orientations of the spins, but we should 
make a more concrete connection between this 
idea and the original presentation of the model. 
We would like to remind the reader that even 
though the angles Oi of the pendula are small, 
the long suspensions of the bobs (£) allow fa to 
cover the entire system which, rescaled, has di¬ 
mensionless length 27 t. The system is also peri¬ 
odic, so the bobs can be thought of as moving 
on a unit circle where the position of the i th bob 
is Xi = 2m /N + fa. In order to think of fa 
as the spin orientations, we start by considering 
each bob as living on its own individual unit cir¬ 
cle. An example of these unit circles is shown 
in Fig. [2j a visual aid to the following. Imagine 
stacking horizontal circles in the vertical direc¬ 
tion and rotating each by an angle 2k/N from 
the one below. The projection of these circles 
onto the horizontal plane would be the system 
viewed in x, i.e. the HMF model. If we twist 
the stack so there is no rotation between adja¬ 
cent circles and then project onto the horizontal 
plane it creates the picture viewed in </>, where 
the pivot points are all aligned. The reason for 
this artificial construction of stacked circles is 
partly to pictorially depict the transformation 
between x and 0 and partly to show how rhi (as 
defined) is just the orientation of the z th spin 
in the <f> picture. Said differently, each circle in 
the cj) picture represents a spin with an orienta¬ 
tion in the horizontal plane determined by fa ; 
an infinite-range classical mean held spin model 


described by Eq. (15). 


j th spin pair depend on both their relative ori¬ 
entation as well as the difference between their 
indices. In the following discussion it will some¬ 
times be convenient to speak about 1 fa in the spin 
language. 

We will prove that in the 4> picture, the sys¬ 
tem in equilibrium with a heat bath is equiva¬ 
lent to the HMF model (the x picture) in equi¬ 


III. EQUILIBRIUM 

In this section, we solve for the canonical par¬ 
tition function, in the (j) coordinate frame using 


the Hamiltonian in Eq. (15) and show that, in 


equilibrium, the HMF model describes the an- 
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a) X b) 0 



Figure 2. Example of a system of N = 8 particles 
when viewed as individual spins in the a) x coordi¬ 
nate frame, and b) the <j> coordinate frame, a) In the 
x coordinate frame the direction of the i th spin given 
by the angle Xi is expressed as Xi = 2iri/N + fa. Al¬ 
ternatively, Xi could be thought of as the position of 
* th particle on the unit circle, shown at the bottom 
of the figure as the projection of all positions onto 
the horizontal plane. The black circles on the rings 
in the figure mark the location of the pendulum piv¬ 
ots at 2ni/N in x. b) Twisting the column of rings 
in a) such that the pivots are aligned transforms the 
system into the fa coordinate frame. In this picture, 
the direction of the i th spin is just given by the angle 

fa- 


gles of long pendula with long-range interacting 
bobs. In order to solve the configurational piece 
of the partition function the Hamiltonian must 
be modified. Using the cosine and sine sum and 
difference identities twice, the potential interac¬ 
tion piece of the Hamiltonian Hj can be written 
as 




cos 


h3 


27 t ( i-j) 


[cos fa COS 4>j 


-sin 


2n(i-j) 

N 


N 

+ sin fa sin fa] 

[sin fa cos fa — cos fa sin fa] 


(16) 


The coefficients in the Hamiltonian 
cos [27 x(i — j)/N] and sin [2w(i — j)/N] should 
be thought of as matrices with components Cij 
and Sij respectively. The Hamiltonian can now 
be written in the form 


T/\r 2J cos cos fa + sin faCij sin 4>j 
ij 

— sin fa cos <f>j + cos fa Sij sin fa), (17) 

which is suggestive because it can be regarded 
as the matrix equation 


rr _ J_ 

1 2N 


(cos fa, COS fa, ..., COS (/>n)C 


( cos fa \ 
cos fa 


+ (sin fa, sin fa, ...,sin 4>n)C 


— (sin fa, sin fa, ...,sin (/)n)S 


+ (cos fa, COS fa, ..., COS 4>n)S 


\COS (f> N J 

( sin fa \ 

sin fa 

\sin (j) N J 

( cos fa \ 
cos fa 

\cos 4>n) 

( sin fa \ 
sin fa 

\sin (f> N J 

(18) 


It is helpful to consider the particles positions on 
the unit circle with respect to their pivot (fa as 
magnetizations. Defining 
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rhi = (cos <j>i, sin </>*), (19) 

and with ntf = (m 0|(i , ..., m N -x,n) where 

fj, holds the place of x or y, the Hamiltonian 
becomes 


'Y / 'T' T' 

H r = — (m x Cm x + m y Cm y 

— rriy Sm x + m^Sm y ). (20) 

A closer examination of the structure of the 
coefficient matrices C and S indicates that they 
take the special form of circulant matrices, and 
thus can be simultaneously diagonalized by a 
unitary matrix U. A circulant matrix has the 
form 


a i 

02 

03 

... ajv-i 

Otv \ 

CLN 

Ol 

02 

■ • • a tv—2 

Oat-1 

CLN-1 

OJV 

Ol 

■ • • a tv—3 

Oat-2 

03 

CI4 

05 

... ai 

02 

02 

03 

04 

... Oat 

Ol / 


( 21 ) 


a special kind of Toeplitz matrix, where each 
subsequent row is a cyclic permutation of the 
row above or below it. Any matrix A with ele¬ 
ments a,ij that can be written in terms of some 
function f(i — j ) is a circulant matrix. Because 
a circulant matrix is a normal matrix it can be 
diagonalized by a unitary matrix. We show that 
C and S are simultaneously diagonalizable by 
showing that they commute, i.e. [C, S] = 0 
where [C. 5] = CS — SC. Starting with the 
second term, —SC = S T C T = ( CS) T which 
is found by arguing that C is symmetric since 
cosine is an even function and does not change 
under the exchange of i and j, whereas S is odd 
because sine is an odd function and does change 
sign under exchange of i and j. The commuta¬ 
tion becomes [C. S] = CS + (CS) T . Also, an 
odd function multiplied by an even function re¬ 
sults in an odd function so the entire matrix CS 
is odd. Therefore ( CS) T = — CS bringing us to 
the final expression [C, 5] = CS — CS = 0. We 
have shown that C and S can be simultaneously 
diagonalized by U. The matrix U is known for 
circulant matrices and is called the Fourier Ma¬ 


trix. 

The matrices C and S can be rewritten as 
C = WD C U and S = U^D S U, where D c and 
D s are diagonal matrices with diagonal elements 
that are the eigenvalues of C and S, respectively, 
which we denote as Xf and Xf. From here on we 
label the indices i from 0 to N — 1. It is worth 
pointing out that due to S being antisymmetric, 
U must be complex. Equation (22) becomes 


Hr = ^(mF x U^D c Um x Am T y U^D c Um y 

- m T y U ] D s Um x + mlU ] D s Um y ) . (22) 

We will move back to the index notation using 
the following relations: 

D c ’ s = A f’%, (23) 


N 

= ( 24 ) 

k=1 

and 

N 

(25) 

k=1 

where X is not to be confused with x. Using the 
Kronecker delta, we set all i = j since these are 
the only nonzero terms. The Hamiltonian is now 
given by 


N—l 


3=0 


(rai 2 A?+iwii 2 a 


+ X*r,Af), ( 26 ) 


with 11XII = XX* and ||y|| = YY*. The in¬ 
clusion of the eigenvalues X c and X s simplifies 
the Hamiltonian further. We will now solve for 
X c and A 5 . Looking at the form of a circulant 


matrix shown in Eq. (21) reminds us that the 


elements of a circulant matrix can be defined 
with a single label. We write the single labeled 
elements of the cosine and sine matrices respec¬ 
tively as 


C; = cos 


2irl 

W' 


(27) 
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and 

2irl 

s i = sm—, (28) 

where l = 0,1,21. The eigenvalues, 
A a , of a N x N circulant matrix A can be 
written in terms of the single label elements 
ai . The j th eigenvalue of A is known to be 
ai exp (2mlj/N), where i is T 
(not an index) and Z = 0,1, 2,..., N — 1. There¬ 
fore, 


*F = E“*(^V Mi/K . < 29 > 

i=o v ' 

and 

A f = E sin (^)e 2 '« /,v . (30) 

2=0 v 7 

Writing cosine and sine in their exponential 
forms gives 


\C 


1 

2 


N-l 

^ |^ e *27rKl+ 1 )/ Ar 
2=0 


_)_ g*27r/(j—1)/Af 


, (31) 


. Af-l 

Aj = Z! ^ J e *27r/(i+l)/iV + e i2nl(j-l)/N 
2 2=0 

(32) 

The above representations of the eigenvalues 
show that C and S each have only two non-zero 
eigenvalues corresponding to j = 1, N — 1 given 
by Af = A£_ x = IV/2 and Af = (A®_ x )* = 
iN/2. The Hamiltonian simplifies greatly to 


H : = l- {WXx + iTi || 2 + ll-XTjv-r - il)v-i|| 2 ) • 

(33) 


The representation of Hi in Eq. (33) must 


be further modified before the partition function 
can be found. We do this by splitting the Fourier 
matrix U into its real and imaginary compo¬ 
nents, ciik and bik, given by 


1 


a ik 


<*> 


: COS 


and 


1 (2-rrik 

bik = —^= sin —— 

Vn \ N 


(35) 


This was done to write the absolute squares 


in Eq. (33) in terms of the squares of di k and 
bik■ By noticing that a lk = a(jv-i)fc and b lk = 
—b(N-i)k we write the configurational partition 
function as 


Zj = A j d N (j)e^&*i aikm Z~ blkrn k}y 

x e ^(E k[ b lkml+a lk m y k }) 2 ^ 


where f5 = 1/ksT. 

The Hubbard-Stratonovich transformation is 
now applied twice, once to each quadratic quan¬ 
tity in the partition function. The integration 
variables introduced through this transforma¬ 
tion are z± and z 2 with subscripts for first and 
second quadratic quantities, respectively. After 
after switching the order of integration, we find 


Z T = 


A 


27 r/Ty J_ c 


dzicfeie-tf+^ )/2 * J! 


d(j) k c 


(zxa lk +Z2b lk )cos(j) h +{z2a lh —zib lk )sin <j> k 


(37) 


The integration can be performed using the iden¬ 
tity 



dcjie^ cos ^ +v sin ^ 


2vr/ 0 (v / e 2 + r/ 2 ) (38) 


where 

f+V 2 = ( z 1 d lk +z 2 bi k ) 2 + (z 2 d lk -zibi k ) 2 (39) 
which simplifies when a and b are included to 


zi 


Vn 


cos 


2nk 

~w 


+ Z 2 


Vn 


sin 


1 

f 2itk\ 

Z 2 . - cos 1 

L Vn 

\ N ) 


- zi 


Vn 


sm 


2irk 

~w 

2irk 

~W 


= -jjVl + zl) 


(40) 
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It is convenient to make a change to polar coor¬ 
dinates by introducing z = \f z\ + z \, following 
which the partition function can be written as 




Vn 


(41) 


Equation (41) is recognized to be an intermedi¬ 


ate step of the solution to the canonical parti¬ 
tion function for the HMF model. From here we 
jump to the main results, the details of which 
are included in the HMF literature mum ■ 


The integration over z in Eq. (41) can be 


preformed using the saddle point approximation. 
The rescaled free energy per particle follows as 


PF = -?- + inf 

Z z 


—z 


2/3 


ln27r Io(z) 


(42) 


The expression above permits a convenient path 
to finding the phase transition. Solving for the 
minimum values of z in order to satisfy the last 


term in Eq. (42) results in the equation 


z 

P 


h{z) 

h(z) 


= 0 , 


(43) 


which can be solved self consistently for 2 and 
represented graphically for different values of (3 
as in Fig. |3j The reader will see that after (3 is 
increased passed the critical value ((3 = 2) there 
are two well-defined solutions. 

The Hubbard-Stratonovich transformation 
decouples spin-spin (squared terms in the Hamil¬ 
tonian) contributions to the partition function 
at the price of needing to create a linear inter¬ 
action between each spin with an auxiliary held 
z\n\. Again, a more detailed procedure can be 
found in [6, .91 [15] where discussion of the inter¬ 
nal energy in the equilibrium state is followed 
by non-equilibrium behavior of the system pre¬ 
pared in microcanonical ensembles. Here we will 
simply touch on the most important point of the 
equilibrium behavior, being that for (3 < 2 the 
system is paramagnetic but for (3 > 2 a pitch- 
fork bifurcation occurs resulting in two stable 
solutions. At this point there is a discontinuity 
in the free energy, a second order phase transi¬ 
tion occurs and the system can maintain finite 
magnetization. In this case, the order parame- 



Figure 3. The solid (black) curve is the fraction of 
modified Bessel functions I\(z)/Ifaz), dashed (green) 
is z/(3 for (3=1, solid (red) line is z/(3 for beta = 2, 
dotted (blue) is zj(3 for (3 = 4, all as a function of z. 
(blue) line is z/(3 for (3 = 4. The values (3 = 1,2,4 
correspond to the pre-phase transition, critical, and 
post-phase transition values in that order. 


ter is the total magnetization M = 
where rhi was defined to be (cos fa, sin fa). 

Showing that the canonical partition func¬ 
tion in the 0 coordinate frame model and the 
HMF model are equivalent necessitates a more 
detailed discussion of the equilibrum behaviors 
in the </> frame. Carnpa et al. |9j, i R their re¬ 
view of the HMF model, rigorously show ensem¬ 
ble equivalence between the canonical and mi- 
crocanoical ensemble of the HMF model. In light 
of this fact, a large N microcanonical simulation 
should be able to produce equilibrium behavior 
like the phase transition mentioned above. The 
temperature in a numerical simulation of a sys¬ 
tem with many particles would be “set” through 
a choice of the initial momenta distribution. In 
this type of simulation, it is common practice 
to compute the order parameter and free energy 
HIM, begging the question: does a large mi- 
crocaonical simulation of Eq. (15) approximate 
the expected equilibrium behavior? Also, since 
the index-dependent model in equilibrium with 
a heat bath can be described by the HMF model, 
would the dynamics of such a simulation quali¬ 
tatively resemble those in the HMF model? The 
answer to both of these questions is no if one 
were to find the equations of motions in fa for 
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some large N and then compare them to an HMF 
model or the x coordinate frame. As stated, this 
discrepancy may appear to detract from our re¬ 
sult. Indeed, it uncovers a conceptual omission 
in the model, but it is one whose rectification 
gives insight into the models ensemble equiva¬ 
lence property of the model, or lack thereof. The 
omission was in the arbitrary scaling of x which 
we will now rectify. 


We introduce the parameter L which gener¬ 
alizes the scaling in Eq. ([ 2 ]) to 


2itL 

~Nd“ 


making the position of the i th particle 


2irLi 2irL £ 

- 77 T AT 7 ') ■ 


N 


N d 


(44) 


(45) 


and changing the definition of 4 >i to 4 >i = 

2'KLiOi/Nd. It can be shown that the intro¬ 
duction of L only changes the final result of 
the partition function by a constant factor of L 
due to the enlarged limits of integration. Nu¬ 
merically, we find is that if L » 1, then the 
simulations in cj) closely reproduce the dynam¬ 
ics of HMF model simulations (dynamics in x). 
Therefore, for large L the mirocanonical simula¬ 
tions can approximate equilibrium and the an¬ 
swers to the previous questions - does a large mi- 


crocaonical simulation of Eq. (15) approximate 


the expected equilibrium behavior, and since the 
index-dependent model in equilibrium with a 
heat bath can be described by the HMF model, 
would the dynamics of such a simulation quali¬ 
tatively resemble those in the HMF model? - be¬ 
comes yes. Alternatively, the coordinate frame 
inequivelence is most extreme for small L. These 
numerical results were found using initial con¬ 
ditions that are randomly distributed fa about 
the domain [— Ln, Lit). It should be stated that 
for the rest of this paper we work with L = 1 
becuase we are inetersed in cases where the cj) 
coordinate frame is markedly differnet than the 
x coordinate frame. 


IV. NON-EQUILIBRIUM RESULTS 

For a system of pendula, it is interesting to 
study an initial configuration where all pendula 
are set to random small displacements from fa = 
0. Specifically we initialize the i th pendulum 
angle, fa, randomly in the range [—tt/N,tt/N). 
In x the indices are ordered in x such that 
x\ < X 2 < X 3 < ... < xjv and the i th bob 
is randomly distributed in the range [27rz/IV — 
ir/N, 2iri/N -\-n/N). It is possible to make some 
general statements about the dynamics of this 
configuration in x using the equations of motion. 
Expressing the Hamiltonian with terms that are 
quadratic in cj) yields 


H T = 

TeI™ 

2N ^ 1 

27r(i - 

N 


ij K 



— cos 

’2vr (i-j)' 

1 ^ 
I'-s.bO 

N 


2 



With this expression, the equations of motion for 
the i th particle can be written as 



(47) 


from which we obtain 


ti = —Y 

2 N ^ 


cos 


2 n(i-j) 

N 


i&j fit) 


+ sin 


2 tt(i- j) 


N 


>. (48) 


In the above equation, the last term and the 
cos [27 t(z — j)/N]fa term sum to zero, leading to 


~7 
2 N 


E 


cos 


27r(j - j) 
N 


" 3 ■ 


(49) 


Using the difference formula, we write the accel¬ 
eration as 
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cos I fai) + sin (M2} 

(50) 

where /i i = 4>jCOs(2Tri/N) and /i2 = 
(f)j sin fani/N). The mass (moment of inertia) 
has been set to unity so the above expression 
is the force as a function of index, fa = F{i). 
If fa\) and fa 2) are known, then the initial dy¬ 


namics of the system are elucidated by Eq. (50), 
but in the case of randomly initialized fa the 
(^1) and (^2) are also random and can be dif¬ 
ferent from one another in both magnitude and 
sign. However, a general description of the re¬ 
sults can be given without exactly knowing these 


coefficients. Equation (50) shows that the initial 


force on a given particle depends on its position 
because the indices are ordered in x. In the con¬ 
tinuum (thermodynamic limit), the force takes 
the form 


—T 

F(x) = —(fa 1) cosx + (/r 2 ) sinx). (51) 

Therefore fai) and fa 2) partly play the role of 
the amplitude of this force as a function of x, 
but also can shift the cos x + sin x spatial depen¬ 
dence, which is periodic over the system length. 
In Fig. [4j we show a fit of the force as a function 
of x using Eq. (51) as well as the actual force 


calculated for an example set of initial condi¬ 
tions. The domain in Fig. [4] can be split into 
two pieces (independent of fa)- one where the 
particles experience a positive force, the other in 
which the particles experience a negative force. 
As time is increased, the movements of the parti¬ 
cles evolve the coefficients fa\) and (^2) in such 
a way that the magnitude of the force decreases 
to zero for all particles and then switches sign 
complementary to the original force. This re¬ 
sults in a standing compression wave of the par¬ 
ticles with a wavelength 2tt. The compression 
wave is not stable and eventually two clusters 
form about each node. These two clusters are 
often referred to as a “bicluster”, or the antifer¬ 
romagnetic state in the HMF model, and have 
been explained by Bane et al. by analysing the 
Vlasov equation. They find that the initial com¬ 



Figure 4. Numerical (blue) and theoretical (red) 
value of the t = 0 force felt by each particle as a 
function of its position. This configuration was made 
with N = 200 and 7 = 10. The initial fa used for 
the numerical calculation were chosen randomly in 
the range [0,27 t/N) which was restricted to positive 
values so that the sign of the initial fai) and (/12) 
were known to be positive. This is not significantly 
different than when the range of fa is centred about 
0. The theoretical curve is fitted using Eq. (51) with 
fai) = 0.0109 and fa 2) = 0.0163. 


pression wave (referred to by a different name) 
creates an effective double-well potential giving 
rise to the bicluster m • The question of the 
bicluster stability has not yet been definitely an¬ 
swered, but for a detailed discussion we refer the 
reader to Leyvraz et al. M- Given the simple 
mapping between the <f> and x coordinate frames, 
we should also be able to show the initial form of 
the force in x as well. As presented in Eq. (15), 
the Hamiltonian in the x coordinate frame only 
differs from the HMF model by a constant 7/2. 
In x, Hi is 


N 

Hi = fajyYl cos ( Xi ( 52 ) 

i,j= 1 

Using the difference identity, we fold the equa¬ 
tions of motion for the i th particle to be 

Xi I — sin Xi ^2 cos x j + cos x i ^2 s * n x i 
\ j j 

(53) 

The sums over cosine and sine of Xj play the 
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same role as {/jl i) and (^ 2 ), and the force at a 
given position Xi is clearly of the same form as 
that shown in Eq. (50). 

Depending on (hi) and (//2), all 0* oscillate 
about zero with amplitudes and phases that de¬ 
pend on their location Xi as discussed above. 
As the clustering in x begins, the 0, begin to 
spread out over the full domain [0, 2 n) and con¬ 
tinue to do so until it is covered. The more in¬ 
teresting case in 0 is when all 0j are initially 
randomly distributed in ranges that depend on 
their index, specifically when 0j are chosen in 
the ranges. [2m/N — ir/N,27ri/N + -k/N) so 
that 0i < 02 < 03 < ••• < 0 jv- It should be 
noted that in this new configuration the dynam¬ 
ics in x are nearly identical to the configuration 
previously discussed for ordered Xi . The dynam¬ 
ics in 0 differ drastically between the two cases 
though. In this ordered angle case, we find some 
interesting grouping of the scaled angles. 

Initially the bobs oscillate with an amplitude 
that depends sinusoidally on their position in 0, 
similar to the previous discussion in the x pic¬ 
ture but with four nodes where the 0* remain 
relatively stationary. Once again this behavior 
could be thought of as a standing compression 
wave, but in 0j it has a wave length of ir whereas 
in the x picture it had a wavelength of 27t. As 
the system evolves, all 0j slowly begin to shift to¬ 
wards the nodes of this standing wave until there 
are four clusters of the angles. After some time, 
the angles begin to re-distribute themselves ran¬ 
domly about the domain. The distribution of 0* 
in these three regimes is summarized in three his¬ 
tograms shown in Fig. [5j Aside from the number 
of clusters, there are two primary differences be¬ 
tween the clustering in 0 and the clustering in x: 
(i) The clustering in 0 only occurs when the an¬ 
gles are ordered in the method described above, 
whereas the dynamics in x look identical regard¬ 
less of the distribution in x, presuming it is some¬ 
what homogeneous about the domain, (ii) The 
clustering in 0 is a quasistationary state whereas 
the clustering in x exists for much longer times 
regardless of the system size. Since the cluster¬ 
ing in 0 is quasistationary, a properly prepared 
system could exist in the clustered angle state 
for an arbitrarily long time but only for large 
N. We can view the effect of increasing N and 



0.25, 
! 0.20 
0.15 
0.10 
> 0,05 
O.OO 1 


Figure 5. These three histograms are made by bin¬ 
ning 0j of 50 particles over three different periods 
of time with 7 = 10 . Going from top to bottom 
each period of time belongs to the dynamical regime 
of: standing “compression wave” in 0* from initial 
configuration of 0* £ [i — 2n/N, i + 2tt/N), clustered 
motion about the four initial nodes of the compres¬ 
sion wave, 0j disordered final state. Specifically the 
values of t are: ti = 0, D = 50, t 3 = 100, t 3 = 200, 
t 5 = 7,000, t 6 = 7,200. 


therefore the lifetimes of the clustered states by 
observing the order of the particle index as a 
function of time. In Fig. [6](a)-(c), we show that 
as N is increased, the time it takes for particles 
to fully mix increases. This is shown by plot¬ 
ting the indices on a color scale from 0 (blue) to 
N — 1 (red) along the horizontal axis as time is 
increased along the vertical axis. In Fig. [6](d) , we 
show how the ordering of the particles changes 
at the very beginning of clustering for N = 100. 
Figure. [6] also shows that the compression wave 
is not quasistationary since it quickly reduces to 
the clustered state regardless of N. 


V. CONCLUSION 


Though the application of statistical mechan¬ 
ics and thermodynamics to systems with long- 
range interactions may not always be appropri¬ 
ate, we find that the canonical partition function 
improves our understanding of a system of pen- 
dula with long-range interacting bobs. Solving 
for the canonical partition function of the Hamil¬ 


tonian in Eq. (15), we show that the equilibrium 
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Figure 6. Each particle is colored from blue, rep¬ 
resenting the smallest index, to red, representing the 
largest index. The order of the particles indices at a 
given moment in time is plotted along the horizon¬ 
tal axis. Time increases along the vertical axis.(a) 
N = 15. (b) N = 30. (c) N = 60. (d) N = 100, 
where a smaller range of time is shown in order to 
see the mixing of the indices as the angles begin to 
cluster. 


behavior in the <p coordinate frame is equiva¬ 
lent to the x coordinate frame, i.e. the HMF 
model. As we have argued that the Hamilto¬ 
nian in Eq. ( |15[ ) describes the behavior of the 
angles of repulsive or attracting pendulum bobs 
(see Fig. [Tj) , then the proven equivalence of the 
canonical partition function of Eq. (15) and the 
Hamiltonian mean field model suggests that the 
Hamiltonian mean field model sufficiently de¬ 
scribes the angles of a system of pendula in 
equilibrium. Ensemble equivalence between the 
microcanonical ensemble and the canonical en¬ 
semble is known for the Hamiltonian mean field 


model model [9] and because of this, the micro- 
canonical simulations could be used to approxi¬ 
mate equilibrium behavior. We End numerically 
that in the case of large system lengths, L , the 
dynamics of the system in cp resemble the dy¬ 
namics of the Hamiltonian mean held model, 
equivalently the behavior of the system in x. 
Therefore for large system sizes of long pendula 
in equilibrium, the HMF model describes their 
dynamics and statistics. 

In this paper we also briefly discuss two par¬ 
ticular sets of non-equilibrium results. In one 
case, the system is initialized with small (pi so 
that Xj are distributed relatively evenly through¬ 
out the x domain. This initial configuration es¬ 
sentially gives rise to the “repulsive” low tem¬ 
perature HMF model which exhibits interest¬ 
ing non-equilibrium behavior and is described in 
great detail by naiii]. In the second case, in 
which i pi are ordered by their index i, we show 
there is a compression wave in <p , followed by 
clustering, and finally a mixed index state dis¬ 
playing no apparent order or structure. This is 
in contrast to the dynamics produced by a ran¬ 
domly distributed set of initial (pi which begins 
and then remains in a random disordered state. 
The clustering that can occur in (p is different 
from the clustering in x because it only occurs 
when the angles are initially ordered and because 
it is quasistationary; the lifetime increases with 
the number of particles in the system. 
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